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We present an ab initio framework to calculate anharmonic phonon frequency and phonon lifetime 
that is applicable to severely anharmonic systems. We employ self-consistent phonon (SCPH) theory 
with microscopic anharmonic force constants, which are extracted from density-functional calcula¬ 
tions using the least absolute shrinkage and selection operator technique. We apply the method 
to the high-temperature phase of SrTiOs and obtain well-defined phonon quasiparticles that are 
free from imaginary frequencies. Here we show that the anharmonic phonon frequency of the an- 
tiferrodistortive mode depends significantly on the system size near the critical temperature of the 
cubic-to-tetragonal phase transition. By applying perturbation theory to the SCPH result, phonon 
lifetimes are calculated for cubic SrTiOs, which are then employed to predict lattice thermal con¬ 
ductivity using the Boltzmann transport equation within the relaxation-time approximation. The 
presented methodology is efficient and accurate, paving the way toward a reliable description of 
thermodynamic, dynamic, and transport properties of systems with severe anharmonicity, including 
thermoelectric, ferroelectric, and superconducting materials. 
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I. INTRODUCTION 

Lattice anharmonicity plays an important role in 
characterizing various physical properties of solids and 
molecules, including the temperature-dependence of vi¬ 
brational frequencies, thermal expansion and phase sta¬ 
bility of solids [1]. It is also responsible for the finite 
phonon linewidth and the lattice thermal conductivity 
kl, which is a key quantity when optimizing the thermo¬ 
electric figure-of-merit ZT [2]. The magnitude of anhar¬ 
monicity varies significantly for different materials. For 
example, covalent materials such as silicon, diamond and 
graphene are very harmonic and show high thermal con¬ 
ductivities [3, 4]. Conversely, thermoelectric and ferro¬ 
electric (FE) materials often show severe anharmonicity, 
demonstrated by inelastic neutron scattering spectra and 
ultralow kl values [5-7]. Anharmonic effects can also be 
significant in superconductors [8-10] and materials under 
extreme conditions [11, 12]. To develop a robust under¬ 
standing of anharmonic properties of solids, a reliable and 
versatile computational method is required. Therefore, 
the development of first-principles methods to calculate 
anharmonic properties of solids and molecules has been 
the subject of intense research in recent years. 

Many-body perturbation theory is one approach for 
treating lattice anharmonicity. This technique consid¬ 
ers the anharmonic effects as self-energies [13]. The self¬ 
energies can be calculated using a systematic approxima¬ 
tion to the Feynman diagrams, where the lowest-order 
approximation is usually employed in the ab initio calcu¬ 
lations based on density-functional theory (DFT). Per¬ 
forming this calculation requires the cubic and quartic 
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force constants, which are the third- and fourth-order 
derivatives of the Born-Oppenheimer potential energy 
surface, respectively. The third-order terms can be ob¬ 
tained efficiently and systematically using either den¬ 
sity functional perturbation theory (DFPT) [14] or the 
finite-displacement approach [15]. Using the cubic terms, 
phonon linewidth can be obtained by evaluating the bub¬ 
ble diagram [Fig. 1(b)]. This type of calculation has been 
performed to predict the lattice thermal conductivity of 
many solids [3, 4, 16, 17] and can also be applied to com¬ 
plex materials [18]. To estimate the phonon frequency 
shift due to lattice anharmonicity, one also needs to com¬ 
pute the loop diagram [Fig. 1(a)] using the quartic terms. 
The calculation of the quartic terms can, in principle, be 
achieved using the finite-displacement approach. How¬ 
ever, since the number of quartic parameters increases 
rapidly as the number of atoms in the supercell increases, 
such calculations have only been reported for simple sys¬ 
tems [19, 20]. 

The perturbative approach is valid only when the an¬ 
harmonic self-energies are sufficiently small compared 
with the harmonic frequency. Therefore, one cannot ex¬ 
pect this technique to yield accurate results for severely 
anharmonic systems. High-temperature phases of FE 
material are typical cases where the perturbation ap¬ 
proach fails because of the imaginary frequencies of har¬ 
monic phonons. To overcome this limitation, it is nec¬ 
essary to employ a non-perturbative approach to treat 
anharmonic effects. 

Methods based on ab initio molecular dynam¬ 
ics (AIMD) can consider anharmonic effects non- 
perturbatively. From the velocity-velocity autocorrela¬ 
tion function calculated using the trajectory of an AIMD 
simulation, one can obtain the vibrational density of 
states with full anharmonicity. To obtain the anhar- 
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monic frequency and linewidth of individual phonons, the 
velocity should be projected onto the phonon eigenvec¬ 
tor [12]. Inherent in this procedure is the assumption 
that the phonon eigenvectors are not altered by anhar- 
monic effects. Such an assumption, however, is valid 
only for simple systems containing a few atoms in the 
primitive cell. The temperature-dependent effective po¬ 
tential (TDEP) method [21] is another AIMD-based ap¬ 
proach. The TDEP method optimizes the effective har¬ 
monic force constants within an AIMD simulation at a 
target temperature. This method should be useful in high 
temperature because it allows both the phonon eigenvec¬ 
tors and the internal coordinate system to be changed by 
anharmonic effects. However, since the AIMD is based on 
the Newton equation of motion, the MD-based methods 
cannot account for the zero-point vibration. Therefore, 
these methods cannot be applied to superconductors and 
ferroelectric materials in the low-temperature range. 

Self-consistent phonon (SCPH) theory [22] is another 
approach for including anharmonic effects beyond per¬ 
turbation theory that considers the quantum effect of 
phonons. Other first-principles methods are able to 
compute anharmonic phonon frequencies related to the 
SCPH theory: self-consistent ab initio lattice dynamics 
(SCAILD) [23] and stochastic self-consistent harmonic 
approximation (SSCHA) [24]. To avoid the cumbersome 
calculation of quartic force constants, these methods em¬ 
ploy real-space stochastic approaches and displace atoms 
in the supercell to model anharmonic effects. 

In this study, we have developed an efficient first- 
principles method to treat lattice anharmonicity. The 
method is based on the SCPH theory, and the anhar¬ 
monic frequency is estimated from the pole of the Green’s 
function. The cubic and quartic force constants necessary 
for the present SCPH calculations are extracted from the 
DFT calculations using the recently proposed compres¬ 
sive sensing approach [25]. By combining the perturba¬ 
tion theory with the solution to the SCPH equation, we 
can also estimate the phonon lifetime and lattice thermal 
conductivity of severely anharmonic materials. 

To confirm the validity of our approach, the method is 
applied to the high-temperature phase of SrTiOa with 
cubic symmetry (c-STO). SrTiOa is one of the most 
studied perovskite oxides and is known to undergo the 
cubic-to-tetragonal phase transition at 105 K accompa¬ 
nied by the freezing-out of the antiferrodistortive (AFD) 
soft mode [26-28]. The EE phase transition is not ob¬ 
served, even at 0 K, because of the zero-point vibration. 
Our approach can describe the temperature dependence 
of the soft-mode frequencies and lattice thermal conduc¬ 
tivity of the severely anharmonic c-STO. 

This paper is organized as follows. First, we introduce 
the SCPH theory and details of our implementation in 
Sec. H. We describe the details of the computational con¬ 
ditions, including the compressive sensing of force con¬ 
stants in Sec. HI. The method is applied to c-STO and the 
results are presented in Sec. IV. In Sec. IV B, we exam¬ 
ine the size- and temperature-dependence of anharmonic 


phonon frequencies and compare these results with ex¬ 
perimental values. We also calculate the lattice thermal 
conductivity of c-STO in Sec. IV C to show the validity of 
our approach. Finally, we conclude this work in Sec. V. 

II. SELF-CONSISTENT PHONON THEORY 
A. Potential energy expansion 

The dynamics of interacting ions within the Born- 
Oppenheimer approximation are described by the Hamil¬ 
tonian H = T U, where T is the kinetic energy and U 
is the potential energy of the system. When U is an an¬ 
alytic function of atomic displacements from equilibrium 
positions {u} , it can be expanded as a Taylor series with 
respect to u as 

[/ = C/o + C /2 + C/3 + C/4 + • • • , (1) 

C/n — T ^ ^• • ■ 5 ^n^n) 

n\ 

{£,«:,//} 

^ ( 2 ) 

Here, Un is the nth-order contribution to the po- 
tential energy, uffin) is the atomic displacement of 
the atom k in the £th cell along the ^ direction, 
and (CiKi;... ffn^n) is the nth-order interatomic 

force constant (IFC). In Eq. (1) the linear term Ui is 
omitted because atomic forces are zero in equilibrium. 

In the harmonic approximation, only the quadratic 
term U 2 is considered and cubic, quartic, and higher- 
order terms are neglected. This allows the Hamiltonian 
Hq = T -\-U 2 to be represented in terms of the harmonic 
phonon frequency w. To compute the phonon frequency 
w, one needs to construct the dynamical matrix 

D,ffnn’-q) = ^ ( 3 ) 

where is the mass of atom k, £'k') are the 

harmonic IFCs, and r(£) is a translation vector of the 
primitive lattice. By diagonalizing the dynamical matrix, 
one obtains harmonic phonon frequencies as 

( 4 ) 

where the index j labels the phonon modes for each crys¬ 
tal momentum vector q and Sqj is the polarization vector 
of the phonon mode qj. 

B. Dyson equation 

To derive the SCPH equation, we employ the 
many-body Green’s function theory. The one-phonon 
imaginary-time Green’s function is given as 

^qj,qi'i'^) — ^ 

= Z-^TT{e-f^^Tr[AqffT)Al,{0)]}, (5) 
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where Tr is the time-ordering operator, Aqj(T) = 
is the displacement operator in the 
Heisenberg picture, Z = Tre“^^ is the partition func¬ 
tion, and /? = 1/fcT, where k is the Boltzmann constant 
and T is the temperature. The displacement operator 
is defined as Aqj = bqj + b^^qj where bqj and bqj are 
the annihilation and creation operators of the phonon 
qj, respectively. It is straightforward to show that 
the Green’s function satisfies Gqjj' (r) = Gqjj' (r -I- /3h) 
for — /?/i < T < 0 and Gqjy(T) = Gqjj/ij — fih) for 
0 < r < ph, where we simply denote Gqj^qjr as Gqjj/. 
Because of these properties, we can also show the fol¬ 
lowing result for the Fourier transform of the Matsubara 
Green’s function; 

pPh 

Gqjf{iUJ^)= dTGqjj'{T)^‘^^^, (6) 

Jo 

where u)m = 21 : 11 X 1 /ph is the Matsubara frequency. To 
obtain the Green’s function for anharmonic systems, we 
need to solve the Dyson equation. When one obtains 
Gqjji (iujm) within some approximations, it is possible to 
obtain the retarded Green’s function Gqjjr (oj) by analytic 
continuation to the real axis as Gqjj' (aj) = Gqjj/ (iuJm 
uj + ie) with a positive infinitesimal e. The function Gqjf 
has a pole at the energy corresponding to the renormal¬ 
ized frequency Dqj. In the case of the harmonic ap¬ 
proximation, one can readily obtain the expression for 
^qjj' (^) 




2u)n 


I _ X.., 

,2 

%3 


(7) 


Therefore, the free-phonon Green’s function is diagonal 
in the phonon polarization index j and can be obtained 
from the harmonic phonon frequencies. 

To estimate the phonon Green’s function Gqjj/ (oj), and 
thereby obtain the anharmonic frequency flqj, we solve 
the Dyson equation 

[GqM]-i = [G°qM]-i-Sq(a;). (8) 

Here we denote the retarded Green’s functions in the ma¬ 
trix form and Sq(uj) is the phonon self-energy, which can 
be estimated within a systematic diagrammatic approxi¬ 
mation. Since the left-hand side of Eq. ( 8 ) becomes zero 
at the frequencies of the renormalized phonons, finding 
the solution {Dq^} is equivalent to solving the following 
equation 


det{[G°H]-i-Sq(a;)} = 0. (9) 

By multiplying det (Aq ) from the left and right of Eq. (9) 
with the diagonal matrix Aqjji = 2ujqjSjj', one obtains 
the following SCPH equation: 


This equation needs to be solved self-consistently because 
the self-energy is a function of the solution uj. In the 
present study, however, the w-dependency in Eq. (11) 
can be neglected because we consider only the first-order 
contribution to the phonon self-energy Sq'*^, which is in¬ 
dependent of UJ, as will be described in Sec. HG. Nev¬ 
ertheless, the self-consistency is retained in the SGPH 
approach because the self-energy is a function of phonon 
frequencies and polarization vectors, which themselves 
are updated by diagonalizing the matrix Vq. 


C. Anharmonic self-energy 

Solving the SCPH equation requires a diagrammatic 
approximation to the phonon self-energy Sq(oj). In this 
study, we consider anharmonicity up to the fourth order, 
i.e. H = H 0 + U 3 + U 4 , where C/„ is the nth-order contri¬ 
bution to the potential energy surface expressed in terms 
of the displacement operator A. This can be obtained by 
substituting 


(iK) = (NM,^) /^Hqe^(K; 9 )e 




( 12 ) 


for Eq. (1), where q labels the phonon modes defined as 
9 = (9 j j) and —q = (—q,j), and N is the number of q 
points. We then obtain the following result: 


nl 


A(qi -f • • • -f <7„) 


{9} 

X Aqj • • • Aq^ . 


$(gi;... ;gn) 


(13) 


The function A (< 7 ) becomes 1 if <7 is an integral mul¬ 
tiple of the reciprocal vector G and is 0 otherwise. 
d>(gi;...; ( 7 „) is the reciprocal representation of the nth- 
order IFCs defined by 


d>(gi;... ;(7„) 

= (M„j •••M^„)“5e^^(/ti;(7i)---e^„(K„;g„) 

f 2 ,...,G 

(14) 


In solving the SCPH equation, we consider only the 
first-order contribution to the phonon self-energy due to 
the quartic term 



1 fi^iqj-,-qj'-,qi-,-qi) 

2 4,^UJqjUJqjiUJq^ 

x[l + 2n(wqj], (15) 


det — \^(a;)} = 0 , ( 10 ) 

^qjj' ~ ~ {‘^^qj)^ i/^^qj') ^ ^qjj' (^)- (H) 


which corresponds to the loop diagram shown in 
Fig. 1(a). Here, n(oj) = — 1]“^ is the Bose-Einstein 

distribution function. Since we continue the iteration 
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(a) qiji 


(b) qiji 



FIG. 1. Diagrams of the self-energies considered in this study. 

(a) The first-order diagram associated with the quartic term. 

(b) The second-order diagram associated with the cubic term. 


cycle of the self-consistent equation [Eq. (11)] until we 
obtain a convergence with respect to the anharmonic fre¬ 
quencies, the SCPH equation automatically includes an 
infinite class of anharmonic self-energies that can be gen¬ 
erated from the loop diagram. In this study, we consider 
the off-diagonal components of the self-energy to allow 
for polarization mixing (PM), which we found to be im¬ 
portant for c-STO, as will be discussed in Sec. IV B. If 
we neglect the off-diagonal elements, Ri , 

the SCPH equation can be reduced to the diagonal form 


= uj^ 
q q 




j{o) 

Q 



M(g;-g;qi;-qi) 


[I + 2n(f2qj)]. 


(16) 

(17) 


This equation is equivalent to the one derived by a varia¬ 
tional approach where the anharmonic free-energy within 
the first-cumulant expansion is minimized with respect to 
trial frequencies [29]. 

To calculate the phonon linewidth, one needs to con¬ 
sider the bubble self-energy shown in Fig. 1(b), which is 
the contribution from cubic anharmonicity given as 

yib) N _ J_ ^4»(-gj, qi,q 2 )^{qj', -gi, -qi) 

..T - 2N 

X A(-g + qi-bq2)-F(*u;m,l,2). (18) 


Here we introduced the oj-dependent function T defined 
as 


J'(iwm,l, 2 ) = ^ 


cr=±l 


1 + ni +n2 
iujm + tT(a;i -I- W 2 ) 

^ n2-ni 

iujm + cr(a;i - a;2)J ’ 


(19) 


where we symbolically denote n{u!q^) and Wg. as rii and 
uji, respectively. We will consider the contribution from 
this diagram in a perturbative manner whereby the equa¬ 
tion (18) is evaluated using the phonon frequencies and 
polarization vectors obtained as a solution to the SCPH 
equation. It should be noted that there is another second- 
order diagram that contains two four-phonon vertexes. 
Although we do not consider that contribution for com¬ 
putational reasons, it is, in principle, possible to extend 
the theory to include higher-order diagrams [30]. 


D. Computational implementation 


In this section we describe the details of the compu¬ 
tational implementation used to solve the SCPH equa¬ 
tion efficiently. The most expensive part of the SCPH 
equation is the calculation of the quartic coefficients in 
Eq. (15), which are changed in each cycle of the itera¬ 
tive algorithm through an update of the phonon eigen¬ 
vectors. To avoid recalculating the quartic coefficient in 
each cycle, we employ a unitary transformation of the 
eigenvectors, as will be described below. Our approach 
is inspired by the method proposed by Hermes and Hi- 
rata for molecules [31], which we extended to periodic 
systems at finite temperatures. 

First, we construct the dynamical matrix D{q) from 
the harmonic IFCs and calculate eigenvalues and eigen¬ 
vectors {ujq, g)} for the gamma-centered W x ^^2 x 
iVa g-point grid. We then calculate the matrix elements 
Fqqi,ijke = ^'(gf;-gj;gifc;-gi£) by Eq. (14) using the 
harmonic eigenvectors and quartic IFCs. Here, the index 
g is restricted to the irreducible points that are com¬ 
mensurate with the supercell size, whereas the index gi 
includes all of the iVi x A 2 x A 3 grid points. The next step 
is to diagonalize the following SCPH equation, which can 
be obtained from Eqs. (11) and (15): 






h[l + 2n{u}q^k)] 

ijkk 


qi,k 


2uj, 


qik 


Here we added the superscript 1 to the matrix Vq to ex¬ 
plicitly show that it is the first iteration of the SCPH 
equation. Then, by diagonalizing the Hermitian matrix 
Vq as Vq = CqWqCl, we obtain the updated phonon 

fll - 

frequencies wT = WTj. The corresponding polariza¬ 
tion vectors can be obtained from the unitary matrix 
Cq. Let Eq and Eq^ denote the s x s matrices de¬ 
fined as Eq = {eqi, ..., Eqs) and eI^^ = ..., el/]), 

respectively, where s is the number of phonon modes. 
It can then be shown that Eq and are unitary 

transformations of each other, which can be written as 
Eq^ = EqCqK Because the phonon polarization vectors 
are modified in this manner, we need to modify Eq. (20) 
for the next iteration of the SCPH equation. The equa¬ 
tion for the nth step of the iteration is given as 

VSl = i ^ Fqq,,,,ulC^;J, (21) 

qi,k,i 


where K. is defined as 

1^, = + (1 - 

„[n] 

a.ki'^a.ki \^ak^ak / 


( 22 ) 


_ ^[n] ^[n]* 

/ V ^q,ki^q,kj 


h[l + 2 n(u;J^y] 


2uj 


qik 


( 23 ) 
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Here we have used the fact that the mean square displace¬ 
ment of normal coordinate Qqj is given as (QqjQqj) = 

2 ^ = 2 l^[l + 2n(wqj)]. In the classical limit 

(/3 —^ 0), the expectation value would be (QqjQqj) = 
kTu}~^. In addition, we introduced the mixing parame¬ 
ter a in Eq. (22) to improve convergence. 

After we obtain Vq'^^ and c]q^ for all irreducible q 
points, we construct the new dynamical matrix as 

r)[n] _ ip[n]-ix^[n] 
q q q q 

= (24) 

where lEq”] = {oJ^q})'^5ij is the diagonal matrix. Using the 
dynamical matrices, we construct dynamical matrices for 
the star of q using the unitary transformation: 

Dfq = r,({5|u(^)})r>Hrt ({^1^(5)}). ( 25 ) 

Here, Tq is the unitary matrix associated with the sym¬ 
metry operation {5|t;(S')} where S is the 3x3 rotation 
matrix and v(S) is the translation vector. The detailed 
expression for can be found in Ref. 32. Finally, we 
construct the dynamical matrix in real space by taking 
the inverse Fourier transformation 

= (26) 

Q 

from which we obtain and for the dense A^i x 
A 2 X N 3 grid points, which are necessary for the next 
iteration of the SCPH equation, by Fourier interpolation. 
For polar semiconductors, the non-analytic part of the 
dynamical matrix is accounted for using the mixed-space 
approach [33]. 

We iterate Eqs. (21)-(26) until convergence is achieved 
for all phonon frequencies at the irreducible q points. 
We initialize the frequency and the unitary matrix as 
= \iOqj\ and = 5ij, respectively. Whenever 

we encounter an imaginary branch, we replace the fre¬ 
quency with its absolute value. After the calculation has 
converged, the anharmonic frequencies and eigenvectors 
for a dense q grid, which are necessary for the subse¬ 
quent calculation of phonon lifetime and lattice thermal 
conductivity, can be obtained by Fourier interpolation. 

III. SIMULATION DETAILS 
A. DFT calculations 

Ah initio DFT calculations were performed using 
the Vienna ah initio simulation package (VASP) [34], 
which employs the projector augmented wave (PAW) 
method [35, 36]. The adapted PAW potentials treat the 
Sr 4s^4p®5s^, Ti 3s^3p®3(i^4s^, and O 25^2^"^ shells as 
valence states. A cutoff energy of 550 eV was employed 
and the Brillouin zone integration was performed with 


the 12x12x12 Monkhorst-Pack fc-point grid. We em¬ 
ployed the PBEsol exchange-correlation functional [37], 
which was reported to work exceedingly well for predict¬ 
ing equilibrium volume and harmonic phonon frequency 
of BaTiOa and SrTiOs [38]. The optimized lattice con¬ 
stant is 3.896 A, which agrees well with the experimental 
value of 3.905 A (Ref. 39, 293 K) and the previous DFT 
result of 3.898 A [38]. The non-analytic part of the dy¬ 
namical matrix is considered in all of the following cal¬ 
culations. We calculated the Born effective charges and 
the dielectric tensor of c-STO using DFPT and obtained 
values of = 6.35, Z*(Sr) = 2.55, 2'*(Ti) = 7.35, 
Z*{0)± = —2.04, and Z'*(0)|| = —5.82, which agree well 
with the previous computational result [40] . Because the 
thermal expansion coefficient of c-STO is very small [39] , 
we neglect thermal expansion effects in this study. 


B. Estimation of force constants 

To compute the harmonic phonon frequency, we ex¬ 
tracted harmonic IFCs using the finite-displacement ap¬ 
proach [15]. The calculation was conducted with a 
2x2x2 cubic supercell containing 40 atoms as in Ref. 12. 
We displaced an atom from its equilibrium position by 
0.01 A and calculated atomic forces for each displaced 
configuration. We then extracted .^'k') by solv¬ 

ing the least-square problem 

$ = argminllA$-Fl]2, (27) 

as implemented in the alamode package [41, 42]. Here, 
$ = [$i, 4>2 ,..., is the parameter vector composed 

of M linearly independent IFCs, E is the vector of atomic 
forces obtained by DFT calculations, and A is the matrix 
composed of the atomic displacements. 

To solve the SCPH equation and estimate the anhar¬ 
monic phonon frequencies of c-STO, one has to prepare 
quartic IFCs. Cubic IFCs are also necessary to esti¬ 
mate phonon linewidth and thermal conductivity, as will 
be discussed in Sec. IV C. In principle, one can extend 
the finite-displacement approach to extract anharmonic 
terms, for which multiple atoms have to be displaced si¬ 
multaneously by an appropriately chosen displacement 
magnitude Art. However, finding an optimal value of Ait 
is not a trivial task, especially when imaginary modes 
exist within the harmonic approximation, as in c-STO. 
We found that the finite-displacement approach with 
Alt = O.l A failed to yield reliable fourth-order IFCs that 
could reproduce the double-well potential of the AFD 
mode. To avoid this issue, one may alternatively em¬ 
ploy the AIMD simulation to sample the displacement- 
force data set. This approach works particularly well for 
simple systems such as Si and Mg 2 Si [41]. However, it 
should be noted that as long as one employs the ordi¬ 
nary least-squares method [Eq. (27)], an overfitting issue 
may arise unless the number of individual reference data 
is fairly large compared with the number of parameters. 
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Recently, Zhou et al. [25] proposed a more robust ap¬ 
proach to estimate anharmonic IFCs. Noting that only 
a small fraction of IFCs has non-negligible contributions 
to atomic forces, they developed the compressive sensing 
lattice dynamics method and obtained the sparse solution 
using the least absolute shrinkage and selection operator 
(LASSO) technique. In the LASSO technique, one solves 
the following equation: 

$ = argmin || A$ — F||| -t- A||$||i, (28) 

where the Li penalty term is added to the least-squares 
equation. Owing to the Li penalty term, one can find 
a sparse representation of the basis function, as demon¬ 
strated by the cluster expansion method and the poten¬ 
tial fitting [43, 44]. In this work, we followed the pro¬ 
cedure of the previous study of Zhou et al. to solve the 
LASSO equation. We initially conducted an AIMD sim¬ 
ulation at 500 K for 2000 steps with the time step of 
2 fs. From the trajectory of the AIMD simulation, we 
then sampled 40 atomic configurations that were equally 
spaced in time. For each configuration, we displaced all 
of the atoms within the supercell by 0.1 A in random 
directions. The atomic forces for the configurations pre¬ 
pared in this manner were calculated using precise DFT 
calculations, from which the matrix A and the vector F 
in Eq. (28) were constructed. The LASSO equation was 
solved using the split Bregman algorithm [43, 45], and 
the optimal value of A was selected from the four-fold 
cross-validation score. To ensure that all of the terms 
in the Li term had the same dimension, we scaled the 
nth-order IFCs and atomic displacement by $ —>■ 
and u — >■ u/uo respectively, with uq = 0.4 oq (~ 0.21 A) 
representing the order of the thermal nuclear motion. 

IV. RESULTS AND DISCUSSION 

A. Anharmonic force constants in cubic SrTiOs 

To find a sparse representation of the basis function 
for c-STO, we first prepared a large parameter vector 
$ that included anharmonic terms up to the sixth or¬ 
der. For harmonic and cubic terms, we included all pos¬ 
sible IFCs present in the 2x2x2 supercell. The quartic 
terms were considered up to third-nearest neighbor shells, 
whereas fifth- and sixth-order IFCs were considered for 
nearest-neighbor pairs. We determined a set of linearly 
independent parameters by considering the space group 
symmetry and the constraints for the translational in¬ 
variance [15, 41]. We hxed the harmonic terms to the 
values determined by the finite-displacement approach 
[Eq. (27)] and employed the LASSO technique for esti¬ 
mating the remaining anharmonic terms. The number 
of linearly independent anharmonic parameters M was 
1053, from which a sparse representation was found by 
Eq. (28). 

Figure 2 shows the magnitude of the anharmonic IFCs 
estimated by solving the LASSO equation. Here, the 



Distance (A) 

FIG. 2. (color online). Absolute values of the third-, fourth-, 
fifth-, and sixth-order anharmonic force constants estimated 
by the LASSO technique plotted as a function of interatomic 
distance. The onsite and two-body terms are indicated by 
circles and the three-body terms are indicated by triangles. 




FIG. 3. (color online). Gomparison of (a) potential energy 
and (b) atomic forces sampled by an individual AIMD sim¬ 
ulation at 300 K. The dashed lines indicate cases where the 
results are identical. 


distance for the IFCs related to more than two atoms is 
defined as the distance of the most distant atomic pairs. 
The absence of onsite force constants for the third- and 
fifth-order IFCs is due to the inversion symmetry of c- 
STO. As shown in Fig. 2, the magnitude of anharmonic 
IFCs decays rapidly with increasing interatomic distance, 
which indicates the locality of anharmonic interactions. 
The terms with the largest magnitude occur at a distance 
of 1.95 A and represent force constants between a Ti atom 
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FIG. 4. (color online). Anharmonic phonon dispersion of c- 
STO at 300 K calculated using the SCPH theory with 8 x 8 x 8 
qi points (solid lines). The dotted lines show the harmonic 
phonon dispersion and the open symbols are experimental 
values at room temperature adapted from Refs. 46 and 47. 


and one of the surrounding O atoms. Among the onsite 
quartic terms, $Tr,Ti,Ti,Ti (m = x,y,z) and 
where v is the direction parallel to the Ti-0 bond, are 
most significant. Compared with these terms, the other 
onsite IFCs, including those of the Sr atom, are one order 
of magnitude smaller. 

The accuracy of the IFCs estimated by the LASSO 
equation was assessed by preparing independent test data 
using an AIMD simulation at 300 K for 2000 steps. 
We then calculated the potential energy [Eq. (1)] and 
atomic forces using the atomic displacements {u} and 
the IFCs {$}. In Fig. 3 we compare the potential en¬ 
ergy U — Uq and the atomic forces obtained from DFT 
and with those calculated from the IFCs estimated by 
LASSO. The model potential well reproduced the DFT 
results for various atomic configurations. The relative er¬ 
rors for the test data were 1.4 and 6.1% for the potential 
energy and the atomic force, respectively, which are as 
small as those reported in Ref. [25]. 


B. SCPH solution 

Using the harmonic and quartic force constants ob¬ 
tained from the finite-displacement and the LASSO tech¬ 
niques, respectively, the SCPH equation [Eqs. (21)-(26)] 
was solved numerically. Since we employed the 2x2x2 
supercell in this study, the q point in Eq. (21) was lim¬ 
ited to the irreducible points on the 2x2x2 grid. We 
changed the qi grid to investigate the convergence of 
the anharmonic phonon frequencies. The mixing param¬ 
eter of a = 0.1 was employed for all temperatures ex¬ 
cept those near the critical temperature of the structural 
phase transition, where a much smaller a was required. 

Figure 4 shows the anharmonic phonon dispersion of 
c-STO at 300 K obtained as the solution for the SCPH 


equation. The phonon frequencies are increased by the 
quartic anharmonicity, evident in the low-energy soft 
modes at the F (0,0,0), R (^,^,^), and M (^,^,0) 
points. We investigated the convergence of the anhar¬ 
monic phonon frequency with respect to the number 
of qi points. The results for the lowest-energy soft modes 
at F, R, and M points are summarized in table 1. Our re¬ 
sults indicate that at least 8 x 8 x 8 qi points are needed to 
obtain convergence and a less dense 2x2x2 qi-point grid 
significantly overestimates the Dq values. This occurs 
because the anharmonic phonon-phonon interaction is 
limited only between the zone-center and zone-boundary 
phonons by the 2x2x2 qi grid. Thus, our numerical re¬ 
sults indicate the importance of including mode coupling 
between longer-wavelength phonons to obtain a reliable 
description of the phonon softening in c-STO. The same 
size-dependence should also be inherent in the real-space 
approaches because the available phonon modes are lim¬ 
ited by the size of the employed supercell. 

We also considered the role of the off-diagonal elements 
of the phonon self-energy that cause PM. In Fig. 5, we 
compare the anharmonic phonon frequencies of two zone- 
center optical modes, labeled TOl and T02, obtained 
using the SCPH equation with and without PM. We 
have shown the SCPH results with 2x2x2 qi points, as 
these results will subsequently be compared with those 
obtained using an MD-based approach. In the SCPH 
equation without PM, we neglect the off-diagonal ele¬ 
ments of the phonon self-energy [Eq. (15)], which is ob¬ 
tained by substituting the unitary matrix Cq in Eqs. (23) 
and (24) with the identity matrix. Therefore, the polar¬ 
ization vectors are fixed to the initial harmonic values. 
Fig. 5 demonstrates that PM is vital to describe the anti¬ 
crossing of the TOl and T02 phonon modes, both of 
which belong to the same irreducible representation F 15 . 
In the case when we neglect PM, an artificial crossing oc¬ 
curs around 500 K and the frequencies significantly devi¬ 
ate from those with PM. Therefore, we conclude that the 
harmonic polarization vectors should not be employed to 
predict anharmonic phonon properties of cubic SrTiOa 
and other perovskite oxides having the same symmetry. 

In Fig. 5 we compare results obtained with the SCPH 


TABLE I. Anharmonic phonon frequency (cm“^) of the soft 
modes at 300 K calculated using the SCPH equation with 
various qi-grid densities. The harmonic phonon frequency is 
also shown for comparison. 


qi points 

Fis (FE) 

R 25 (AFD) 

M 3 

2 x 2 x 2 

144 

69 

103 

4x4x4 

138 

46 

89 

6 x 6 x 6 

136 

39 

86 

8 x 8 x 8 

136 

37 

85 

10 x 10 x 10 

135 

36 

85 

12 x 12 x 12 

135 

35 

85 

Frozen phonon 

58i 

76i 

21 






















Temperature (K) 


FIG. 5. (color online). Temperature-dependence of the an- 
harmonic phonon freqnencies of two F15 modes calculated us¬ 
ing the SCPH equation with and without PM, and with the 
TDEP approach (see the text for details). The TOl mode 
corresponds to the FE mode. The open symbols are experi¬ 
mental values for the TOl mode reported by Servoin et al. [48] 
and Yamada and Shirane [27]. 


method with those obtained with the temperature- 
dependent effective potential (TDEP) method [21]. In 
the TDEP method, atomic displacements and forces are 
sampled by AIMD simulations at a target temperature 
and are then used to extract effective harmonic force 
constants by numerical fitting. In our TDEP simula¬ 
tions, we performed MD simulations using the Taylor 
expansion potential [Eq. (1)] instead of AIMD to re¬ 
duce computational costs. Anharmonic terms up to the 
sixth order were considered and the force constants esti¬ 
mated by the LASSO technique were employed. We con¬ 
ducted the constant-temperature MD simulations with 
the 2x2x2 supercell and the temperature was controlled 
by the Berendsen thermostat [49]. We employed a time 
step of 1 fs and conducted the MD simulations for 50000 
steps at each temperature. The last 40000 steps were 
employed to extract effective harmonic IFCs by least- 
squares fitting [Eq. (27)]. Although the anharmonic fre¬ 
quencies obtained using the TDEP approach are slightly 
smaller than the SCPH results, they agree qualitatively 
with the SCPH results, as shown in Fig. 5. We con¬ 
sider this discrepancy to be reasonable for the following 
two reasons. First, the SCPH results include only anhar¬ 
monic self-energies that can be generated from Fig. 1(a), 
whereas the TDEP includes higher-order anharmonic ef¬ 
fects. Among these higher-order terms, the first-order 
contribution due to the cubic anharmonicity, as depicted 
in Fig. 1(b), should have the largest contribution. We 
found that the effect of the diagram in Fig. 1(b) is to 
reduce the anharmonic frequency for the FE mode. Sec¬ 
ond, since the quantum effect of nuclear motion is not 
considered in the MD simulation, the thermal average 
of the squared normal coordinate (Q*Qq) is underesti- 


FIG. 6. (color online). Temperature-dependence of the 
squared phonon frequency of the R25 mode obtained from 
the SCPH theory with various qi-point densities. The open 
circles are experimental values adapted from Ref. 47. 


mated for temperatures below the Debye temperature in 
the TDEP approach. Therefore, the renormalization of 
anharmonic effects is underestimated in the TDEP ap¬ 
proach, which explains why the deviation from the SCPH 
result becomes larger with decreasing temperature (see 
TOl mode in Fig. 5). 

Figure 6 shows a comparison of the temperature- 
dependence of the squared frequency of the AFD mode 
and experimental measurements [47]. As can be seen in 
the figure, the frequency of the AFD mode is severely 
size-dependent. For the 2x2x2 qi grid, we do not ob¬ 
serve a freezing-out of the AFD mode even at absolute 
zero. When we increase the qi-grid density and allow 
interactions with longer-wavelength phonons, we observe 
the precursor of the freezing-out of the AFD mode at 
temperatures near 200 K. Although the soft-mode fre¬ 
quency does not reach zero in the current simulation with 
finite <71 points, this would occur in the thermodynamic 
limit {N —>• 00) as discussed by Cowley [50]. Above 
approximately 300 K, the temperature-dependence can 
be reliably fitted by the equation D^(T) = a(T — Tc)^. 
Applying this equation to the result obtained using the 
12x12x12 <7i grid, we obtain the Tc of the cubic-to- 
tetragonal phase transition as 220 K. 

For comparison, we have plotted experimental results 
in Figs. 4, 5, and 6 using open symbols. The SCPH equa¬ 
tion reproduces the temperature-dependence of the soft 
modes qualitatively, but not quantitatively, i.e., the fre¬ 
quencies of the FE and AFD modes are overestimated 
and underestimated, respectively. Because the ADF fre¬ 
quency is underestimated, the transition temperature 
predicted is twice as large as the experimental value 
of 105 K. We consider this deviation to be acceptable 
because phonon-related properties of ferroelectric mate¬ 
rials are known to be sensitive to the lattice constant 
and exchange-correlation functional employed [38, 51]. 
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In this study, we employed the PBEsol functional to 
avoid problems inherent to the local-density approxi¬ 
mation (LDA) and the generalized-gradient approxima¬ 
tion with the Perdew-Burke-Ernzerhof parameterization 
(PBE) [52]; LDA tends to underestimate the equilibrium 
volume, whereas PBE tends to overestimate it. How¬ 
ever, our numerical results suggest that PBEsol cannot 
give a quantitative description of c-STO. This issue is ex¬ 
pected to be resolved, at least partially, by employing a 
hybrid functional. Wahl et al. [38] investigated the func¬ 
tional dependence of the harmonic frequency in the EE 
mode of c-STO and reported the results of 29i and 74i for 
the PBEsol semilocal and the Heyd-Scuseria-Ernzerhof 
(HSE) hybrid functionals [53], respectively. Since the 
harmonic frequency changes as < 0, we 

expect that the Fock exchange can increase the depth 
of the double-well potential, thereby decreasing the an- 
harmonic frequency of the EE mode. Wahl et al. also 
reported that the energy gain for the AFD phase was 
smaller in HSE than in the PBEsol functional. This in¬ 
dicates that the depth of the double-well potential for 
the AFD mode can be decreased, and the anharmonic fre¬ 
quency can be increased by using HSE instead of PBEsol. 
Therefore, we believe that the quantitative accuracy of 
the present SCPH results could be improved by employ¬ 
ing a hybrid functional, which will be the topic of future 
work. 

In the present SCPH calculations we have not con¬ 
sidered effects related to the cubic anharmonicity, such 
as thermal expansion, relaxation of internal coordinates 
and intrinsic frequency shifts due to the bubble diagram. 
However, these effects can, in general, become impor¬ 
tant in severely anharmonic systems [50], and should 
be considered, especially when one intends to quanti¬ 
tatively compare theoretical results with experimental 
data. Therefore, extending the present ab initio method 
to include these effects, either perturbatively or self- 
consistently, could be another important direction for 
further research and development. 


C. Lattice thermal conductivity 

The lattice thermal conductivity is a key quantity for 
optimizing the thermoelectric figure-of-merit ZT, and it 
has been the subject of intense theoretical study in recent 
years. To show the validity of our theoretical approach 
based on the SCPH equation, we estimated the lattice 
thermal conductivity of c-STO. For this work, we em¬ 
ploy the Boltzmann transport equation (BTE) within the 
relaxation time approximation (RTA), where the lattice 
thermal conductivity is given as 

<{T) = -^ E C,{TX{TXIT)t,IT). (29) 

Here, V is the unit-cell volume, Cq is the lattice spe¬ 
cific heat, Vq = dftq/dq is the group velocity, and 



Temperature (K) 


FIG. 7. (color online). Temperature-dependence of the lat¬ 
tice thermal conductivity of c-STO. The computational result 
is compared with experimental values reported by Muta et 
al. [54] and Popuri et al. [55]. Lines are shown to guide the 
eye. Inset: Calculated phonon lifetime of c-STO at 300 K. 

Tq = [2r,j(Oq)]“^ is the lifetime of phonon q. The phonon 
linewidth rg(a;) can be obtained from the imaginary part 
of the phonon self-energy that results from the cubic an¬ 
harmonicity [Eq. (18)], which is given explicitly as 

r la;! = — V 

qt q" H H H 

X [{uq' + nqii -I- l)(5(a; — Hg/ — Og//) 

— 2{nq' — nq//)6{uj — flq' + Og")] . (30) 

Here, the matrix element $(g,q',q") is calculated from 
cubic IFCs using Eq. (14) with eigenvectors {e^(K; g)} re¬ 
placed by the solution to the SCPH equation {e^(K; q)}. 
The equations (29) and (30) are identical to those that 
have commonly been employed in the thermal conductiv¬ 
ity calculations except that harmonic phonon frequencies 
and eigenvectors are substituted by anharmonic frequen¬ 
cies and eigenvectors, respectively, obtained using the 
SCPH equation. 

Figure 7 compares the calculated thermal conductivity 
of c-STO with experimental results [54, 55]. The calcu¬ 
lation was conducted using the 8x8x8 qi grid for the 
SCPH equation and the 12x12x12 q grid for the BTE- 
RTA equation [Eq. (29)]. Although we observed devi¬ 
ations in soft-mode frequencies, the calculated thermal 
conductivity agrees well with the experimental results, 
as can be seen in Fig. 7. We expect that the agree¬ 
ment could be improved further by employing a finer q 
grid and using a hybrid functional, although such cal¬ 
culations were not performed because of computational 
limitations. In the Fig. 7 inset, we also show the phonon 
lifetime Tg at 300 K calculated by Eq. (30). The phonon 
lifetimes of c-STO obtained from the perturbation the¬ 
ory [Eq. (30)] are found to be even smaller than those of 
PbTe [16], but the kl value of c-STO is higher. This can 
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be attributed to the large group velocities of phonons, 
especially of TO modes above ~ 100 cm“^, which con¬ 
tribute significantly to the total kl value [56]. The life¬ 
time shows a characteristic feature in the low-frequency 
region (< 100 cm“^): the phonon modes split into two 
separate regions in r, > 3 ps and ^ 0.6 ps, where 
the former corresponds to the acoustic modes that fol¬ 
low the frequency dependence of r ^ which has been 
observed in other materials [16, 41], and the latter corre¬ 
sponds to the phonon modes around the R point, which 
indicates the severe anharmonicity of the AFD mode. 

V. CONCLUSIONS 

We developed an ab initio method to compute anhar- 
monic phonon frequencies and lifetimes that can be ap¬ 
plied to severely anharmonic systems. The method em¬ 
ploys anharmonic force constants up to the fourth order, 
which are extracted from DFT calculations using a com¬ 
pressive sensing approach. The frequency renormaliza¬ 
tion associated with the quartic anharmonicity is treated 
non-perturbatively using the SCPH theory. By perform¬ 
ing the perturbation calculation after the SCPH solution, 
we also calculated phonon lifetimes that result from the 
three-phonon scattering processes. 

We applied the method to the high-temperature phase 
of perovskite SrTiOa. Unlike the harmonic phonon dis¬ 
persion, the SCPH solution was free from the imaginary 
branches in the entire Brillouin zone. We found that 
including polarization mixing is important to correctly 
account for the temperature dependence of the phonon 
frequency of the ferroelectric soft mode of perovskite ox¬ 
ides. In addition, we examined the size-dependence of the 
anharmonic frequencies of the soft modes and found that 
long-wavelength phonons significantly reduced the anhar¬ 
monic frequencies, especially for the antiferrodistortive 
mode near the transition temperature. The temperature- 
dependence of the soft mode frequencies calculated using 
the SCPH theory agreed qualitatively well with the ex¬ 
perimental results. However, the quantitative accuracy 
of the present calculations based on the PBEsol func¬ 


tional was unsatisfactory, where we obtained the cubic- 
to-tetragonal transition temperature as Tc = 220 K that 
was twice as large as the experimental value of 105 K. 
Although further theoretical investigations are required 
to understand the origin of this discrepancy, we expect 
that the quantitative agreement can be improved by em¬ 
ploying a hybrid functional. We also calculated the lat¬ 
tice thermal conductivity kl of cubic SrTiOa using the 
Boltzmann transport equation within the relaxation-time 
approximation. The calculated kl values reproduced ex¬ 
perimental results especially in the high temperature re¬ 
gion. The underestimation of kl in the low temperature 
region may be attributed to the overestimation (under¬ 
estimation) of the ferroelectric (antiferrodistortive) soft 
mode, which will be addressed in a future work. 

The present method, which combines the SCPH theory 
with perturbation approach based on anharmonic force 
constants, enables us to obtain the anharmonic phonon 
frequencies and lifetimes at various temperatures effi¬ 
ciently just by changing the occupation number. The 
system size dependency can be investigated using the re¬ 
ciprocal space formalism. Therefore, we believe that the 
present method paves the way for understanding lattice 
anharmonicity and related dynamical and thermodynam¬ 
ical properties of thermoelectric, ferroelectric, and super¬ 
conducting materials. 
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